!*******************************************************************
!> this module has all data and procedures for reading and handling
!  user problem input
!*******************************************************************
module class_input

    implicit none
#include "finclude/petsc.h90"

    !> primary structure for input data
    type, public :: indat
        !> problem input file name
        character(len=80)  :: file
        !> problem description
        character(len=80) :: disc
        !> problem type:
        ! 0 = fixed source (no mult)
        ! 1 = fixed source (mult)
        ! 2 = eigenvalue
        PetscInt :: ptype
        !> number of energy groups
        PetscInt :: numg
        !> number of materials
        PetscInt :: numm
        !> material/group diffusion coefficient
        PetscScalar, dimension(:,:), pointer :: dc
        !> material/group sigR
        PetscScalar, dimension(:,:), pointer :: sr
        !> material/group sigA
        PetscScalar, dimension(:,:), pointer :: ab
        !> material/group nuSigF
        PetscScalar, dimension(:,:), pointer :: ns
        !> material/group xi-sprectrum
        PetscScalar, dimension(:,:), pointer :: xi
        !> material/group sigS
        PetscScalar, dimension(:,:,:), pointer :: sc
        !> material/group sigA
        PetscScalar, dimension(:,:), pointer :: sa
        !> volumetric source by coarse mesh and group
        PetscScalar, dimension(:,:,:), pointer :: src
        !> material placement by coarse mesh
        PetscInt, dimension(:,:), pointer :: mt
        !> number of x and y fine and coarse meshes
        PetscInt :: nxcm, nxfm, nycm, nyfm
        !> x and y coarse mesh
        PetscScalar, dimension(:), pointer :: xcm, ycm
        !> x and y fine mesh (counts)
        PetscInt, dimension(:), pointer  :: xfm, yfm
        !> dx, dy, and dv values
        Vec :: dx, dy, dv
        !> keff iteration precision
        PetscScalar :: epsk
        !> flux iteration precision
        PetscScalar :: epss
        !> maximum number of outer iterations
        PetscScalar :: maxit
        !> left boundary   (0=vacuum, 1=reflect, 2=inc cur)
        PetscInt :: bcl
        !> right boundary  (0=vacuum, 1=reflect, 2=inc cur)
        PetscInt :: bcr
        !> bottom boundary (0=vacuum, 1=reflect, 2=inc cur)
        PetscInt :: bcb
        !> top boundary    (0=vacuum, 1=reflect, 2=inc cur)
        PetscInt :: bct
        !> left incident boundary source
        PetscScalar, dimension(:,:), pointer :: IBSL
        !> right incident boundary source
        PetscScalar, dimension(:,:), pointer :: IBSR
        !> bottome incident boundary source
        PetscScalar, dimension(:,:), pointer :: IBSB
        !> top incident boundary source
        PetscScalar, dimension(:,:), pointer :: IBST
        !> should i print? 1=yes, 0=no
        PetscInt :: printout
        !> where should i print?
        character(len=80)  :: outfile
        !> should i plot? 1=yes
        PetscInt :: plotout ! currently not implemented

    end type indat

    !> unit for output
    integer :: uout

    contains

!************************************************************
!> Subroutine readinput ()
! This subroutine reads the command line input for the file
! containing all the relevant input parameters.
!************************************************************
    subroutine readinput(inp)

        implicit none

        type(indat), intent(out) :: inp
        PetscInt :: io, uinp=5
        ! temporary variables for reading in
        character(80)  :: tmpfile, outfile
        PetscInt :: m,g,gg,numg,numm,nxcm,nxfm,nycm,nyfm,bcl,bcr,bcb,bct,ptype,maxit,printout,plotout
        PetscScalar, dimension(:,:), pointer :: data
        PetscScalar :: epsk, epss
        namelist/matsgeo/numg,numm,nxcm,nxfm,nycm,nyfm,bcl,bcr,bcb,bct,maxit
        namelist/control/maxit,epsk,epss,ptype
        namelist/output/outfile,printout,plotout

        ! get the command input (just the input file name)
        if ( COMMAND_ARGUMENT_COUNT() .lt. 1 ) then
            print *, "*** ERROR: user input file not specified ***"
            return
        else
            call get_command_argument(1,inp%file);
            open (unit = uinp, file = inp%file, action = "read", &
                  status = "old", position = "rewind", iostat = io)
            if (io > 0) stop "*** ERROR: user input file not found ***"
        end if

        ! process input file
        read (uinp, '(a80)') inp%disc
        write(0,*) "initializing problem: "
        write(0,*) trim(inp%disc)
        ! read first namelists
        read (uinp, nml=matsgeo)
        read (uinp, nml=control)
        read (uinp, nml=output)
        ! if we are printing output, open the file
        if (printout .lt. 1) then
            uout = 0 ! easy user option to print to screen instead
        else
            open(uout,iostat=io,file=outfile,status='unknown',form='formatted')
            if (io .ne. 0) stop "*** ERROR: Unable to open output file ***"
        end if
        ! now put these into the input
        inp%numg=numg; inp%numm=numm
        inp%nxcm=nxcm; inp%nycm=nycm
        inp%nxfm=nyfm; inp%nyfm=nyfm
        inp%bcl=bcl;   inp%bcr=bcr;   inp%bcb=bcb;   inp%bct=bct;
        inp%maxit=maxit; inp%epsk=epsk; inp%epss=epss;   inp%ptype=ptype;
        inp%outfile=outfile; inp%printout=printout;  inp%plotout=plotout
        ! allocate the geometry arrays
        allocate( inp%xcm(nxcm+1),         &
                  inp%xfm(nxcm),           &
                  inp%ycm(nycm+1),         &
                  inp%yfm(nycm),           &
                  inp%src(nxcm,nycm,numg), &
                  inp%mt(nxcm,nycm),       &
                  data(numm*numg,4+numg)   )
        ! allocate coefficient arrays
        allocate( inp%dc(numm,numg),       &
                  inp%sr(numm,numg),       &
                  inp%ab(numm,numg),       &
                  inp%ns(numm,numg),       &
                  inp%xi(numm,numg),       &
                  inp%sc(numm,numg,numg)   )
        ! put an allocation error check here
        ! read them

        read (uinp,'(a)')
        read (uinp,*) inp%xcm
        read (uinp,'(a)')
        read (uinp,*) inp%xfm
        read (uinp,'(a)')
        read (uinp,*) inp%ycm
        read (uinp,'(a)')
        read (uinp,*) inp%yfm
        read (uinp,'(a)')
        read (uinp,*) inp%mt
        read (uinp,'(a)')
        if (inp%ptype.ne.1) read (uinp,*) inp%src
        read (uinp,'(a)')
        ! read the cross-section data
        do io = 1, numg*numm
            read (uinp,*) data(io,1:4+numg )
        end do

        ! assign coefficients
        do m = 1, inp%numm
            do g = 1, inp%numg !g+(m-1)*g
                inp%dc(m,g) = data(inp%numg*(m-1)+g,1)  ! diffusion coefficient
                inp%sr(m,g) = data(inp%numg*(m-1)+g,2)  ! removal cross-section
                inp%ns(m,g) = data(inp%numg*(m-1)+g,3)  ! fission cross-section
                inp%xi(m,g) = data(inp%numg*(m-1)+g,4)  ! fission spectrum
                do gg = 1, inp%numg
                    inp%sc(m,g,gg) = data(inp%numg*(m-1)+g,4+gg)
                end do
                if (inp%numg .ge. 1) then
                    inp%ab(m,g) = inp%sr(m,g) - sum(data(inp%numg*(m-1)+1:inp%numg*m, 4+g ))
                else
                    inp%ab(m,g) = inp%sr(m,g)
                end if
            end do
        end do


        deallocate( data )

        !deallocate( inp%dc, inp%sr, inp%ns, inp%xi, inp%sc, inp%ab )
        !deallocate( inp%xcm, inp%xfm, inp%src, inp%mt)


    end subroutine readinput

end module class_input
